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We study a quantum phase transition in a system of dipoles confined in a stack of identical one- 
dimensional ("cigar" -type) layers and polarized perpendicularly to the layers. In this arrangement 
the intra-layer interaction is purely repulsive preventing the system collapse and the inter-layer 
one is attractive. The dipoles may represent polar molecules confined in the optical lattice or 
indirect excitons in multi-layered structures. The transition separates two phases; in one of them 
superfluidity (understood as algebraic decay of the corresponding correlation functions) takes place 
on each individual layer, in the other the order parameter is the product of bosonic operators of 
all layers. We argue that in the presence of finite inter-layer tunneling the transition belongs to 
the universality class of the q = N two-dimensional classical Potts model. For A'' = 2,3,4 the 
corresponding low energy field theory is the model of Zjv parafermions perturbed by the thermal 
operator. Results of Monte Carlo simulations are consistent with these predictions. The detection 
scheme of the chain superfluid is outlined. 

PACS numbers: 



I. INTRODUCTION 

Emergence of Majorana fermions (see in in topological insulators has inspired a search for such 
fermions in other condensed matter systems. In this article we show that parafermions, of which Majorana 
fermions represent a particular case, describe excitation spectra of quantum chains (strings) of polarized 
dipoles. Material realization of such systems has became possible due to the recent breakthroughs in 
creating and trapping high density samples of (polar) molecules 0. As proposed in Ref.Q, multi- 
layered structures of indirect excitons [sf may also form similar systems in the form of excitonic chains. 
Each indirect exciton (not to be confused with the excitons formed at non-F point) has static dipole 
moment due to a spatial separation of electron and holes. Interaction between the dipoles in the A'^- 
layered structure can encourage a formation of the excitonic chains similar to chains of polar molecules. 
Since light field E and excitons are coupled linearly a state of excitonic field ^ is imprinted directly on 
the emitted light. As a consequence, properties of excitonic chains can be explored through light emission 
providing a new powerful experimental tool to study strongly correlated systems. 

So far, quantum chains have been studied in various analytical approximations which neglect tunneling 
of particles along the chains. In Ref.0 it has been proposed that stiff dipolar non-interacting quantum 
chains may form Bosc-Einstcin condensate. Intcr-layerpairing in bilayercd 2D dipolar fcrmionic systems 
has been studied in the ECS approximation in Refs.[7[. The dimerization transition in the 2D multi- 
layered geometry of dipolar fermions was analyzed in Ref.0], and it has also been proposed that for 
strong dipolar interactions long chains can form by the N-clock phase transition Q. Fermionic dipolar 
molecules forming a mixture of single fermions, dimers and fermionic trimers in ID N = 3-layered system 
has been discussed in the ideal gas approximation in Ref.Q. 

In low dimensions quantum fluctuations are enhanced and therefore quantum particles from different 
Id lattices (tubes) may develop strong correlations. This can lead to interesting physics. The difficulty 
is that such system, in general, is not amenable to the standard mean field or perturbation expansion 
methods and one has to resort to a combination of non-perturbative techniques and numerics. A numerical 
study of quantum chains which takes into account the partial-chain exchanges as well as the intra-chain 
dynamics has been performed in Ref.fioj for the case of zero inter-layer tunneling. It has been shown that 
polar molecules in the A^-layered geometry can form flexible (quantum rough) chains, and these chains 
can undergo a quantum phase transition to a superfluid phase characterized by off-diagonal long-range 
order (ODLRO) in the A^-body density matrix, while all M-body density matrices with M < N show 
insulating behavior regardless of the filling factor (provided it is the same in each layer). If the inter- layer 
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FIG. 1: (Color online) Schematic picture of A*' = 3 parallel to each other Id lattices (tubes) stretched along 
X-axis. Each site is depicted by filled circle. Dipole particles are represented by arrows showing polarization of 
their dipole moment (along Z-axis). Chains are understood as bound states of particles from different lattices 
(tubes). The dipoles are allowed to tunnel between nearest sites along X and Z directions. 

(dipolar) interactions arc weak, a stack of N layers features a iV-component superfluid (N-SF). Once 
the interaction becomes stronger, the non-dissipative drag between layers will eventually convert the 
N-SF phase to the flexible chain superfluid (CSF) characterized by ODLRO only in the A^-body density 
matrix [lol |. The corresponding transition is continuous in the ID-geometry and discontinuous in the 
2D-geometry for N >2^. 

In the present work we study the transition between N-SF and CSF states in the presence of finite 
inter-layer tunneling. Since for 2D-layers {N > 2) the transition is of the first order and the inter-layer 
tunneling cannot change this, we concentrate on the ID-geometry depicted in Fig. [1] that is, when layers 
may be considered as tubes. 

Our main findings are the following. A quantum phase transition into the A'^-chain superfluid (in 1-1-1 
dimensions) is in the universality of the classical q = N 2D Potts model. That is, for N = 2,3,4 the 
transition is a continuous one and for > 4 it is of the first order. For = 2 we develop a microscopic 
low energy description in terms of the field theory of two species of Majorana fermions and one gapless 
bosonic field. For A'^ = 3,4 instead of a detailed derivation we present arguments based on symmetry of 
the problem and on results of our numerical calculations. 

The paper is organized as follows. In the Section 2 we discuss possible phases. Then, in Sec. IIIII 
we will describe the microscopic lattice model accounting for a system of N coupled tubes. Then we 
formulate field theoretical description of the system valid in the continuum limit. This description is 
rigorously derived for the case of two tubes (A^ = 2) where the low energy limit leads to a model of 
relativistic Majorana fermions. We also present plausible arguments concerning possible field theory for 
the cases A^ = 3, 4. These arguments are supported by the Monte Carlo calculations presented in Sec. lIVI 
The Monte Carlo procedure is performed for the coarse-grained dual version of the Hamiltonian in the 
discretized time approximation. Finally, in the Conclusion we will give a summary of the main results 
and will briefly discuss the detection scheme. 

II. ORDER PARAMETERS AND UNIVERSALITY OF THE TRANSITION TO THE CSF 

PHASE 

Without inter-layer tunneling in the N-SF state [loj each layer is characterized by its own phase. In 
the case of finite inter-layer tunneling (along Z in Fig. [1]) these phases lock in into a single phase field 
(fi so that the superfiuid (SF) is characterized by the bosonic operator ip ~ e*'^. As is well known, in 
ID the real long range order is substituted by quasi long range order characterized by nonzero stiffness 
and algebraic decay of certain correlation functions at zero temperature. The SF phase of our system is 
characterized by an algebraic decay of the bosonic field. Meanwhile in the Chain Superfluid Phase (CSF) 
correlators of individual Bose operators {x, z)tp{x' , z')) decay exponentially with respect to \x — x'\ 



and an algebraic order pertains to the product of operators of all tubes 
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^jv(x)= n ^ix,z), (1) 

2=l,2,...,Af 

describing the order of quasi-molecular complexes each consisting of N bosons. 

It is important that ^'jv, Eq.(IT]), is invariant with respect to the transformation iplx^z) 
exp{2-Kim{z,x) / N)ip{x, z) where m{z,x) is defined modulo N and obeys the constraint X^z^C-^'^) ~ 
pN,p = 0,1,2,.... Thus, m{z,x) can be broken as m = m' + m into the discrete global part 
m = p, p = 1, 2, A'^ — 1 and the local gauge-type m! {z, x) obeying m' = 0. 

Setting aside the discussion of a possible role of the local-gauge symmetry, we note that the global 
transformation forms a discrete symmetry group which determines the universality of the SF-CSF tran- 
sition. Among the possible candidates one can consider the p-clock model and the (standard) Potts 
model (see in Ref.[Tl[). While the case N = 2 should be assigned to the Ising universality class [l^ . 
the nature of the transition for > 2 is not obvious at all. Naively, one may anticipate the p-clock 
universality because of the nearest neighbor tunneling (between layers). In what follows we will show 
that such expectation is not correct, and the criticality is controlled by the standard Potts model (also 
called as Ashkin- Teller-Potts model). Accordingly, for ID tubes (that is, D = 1 + 1 membranes) it should 
be continuous for iV = 2, 3, 4 and discontinuous for > 4. 



III. MICROSCOPIC HAMILTONAN AND THE EFFECTIVE MODEL FOR BILAYERS IN 

TERMS OF MAJORANA FERMIONS 

Each tube represents an optical lattice occupied by particles with Bose statistics (polar molecules or 
indirect excitons). The microscopic Hamiltonian H describing SF and CSF has the following form: 



x.z 



+ h.c."^ +tj_(^alj^^^^^a^^^ + h.c.^ + ^ X! ^^z;^;'^'"-^^"^;'^' , (2) 



Here aj,^, axz are creation (annihilation) operators creating (destroying) a boson at site x belonging to zth 
layer; Uxz = c-xz^^xz denotes onsite density operator obeying the hard-core constraint; Vxz-x'z' describes 
the matrix element for dipole-dipole interaction between sites (xz) and (x'z'). It is characterized by the 
strength Vd = (ii/b^, where dz stands for the induced dipole moment and bz denotes distance between 
two nearest layers. This interaction is mainly attractive along the z-direction and repulsive along the 
x-direction. 



A. Two chains (A^ — 2). Low energy decsription 

In the low energy limit the microscopic Hamiltonian ([2]) can be replaced by the effective model de- 
scribing the SF to CSF transition in the approximation where the dipole-dipole interaction is reduced to 
the attraction between nearest neighbors layers. Taking into account the single occupancy constraint we 
can rewrite Hamiltonian ^ in terms of the Pauli matrix operators. Restricting ourselves to the simplest 
case of two chains, we have: 

j 2 = 1,2 

+t±Hi''7,2 + h.c.) - Vial^al^ (3) 

where cr^ operators stand, respectively, for the bosonic creation and annihilation ones a^, a in Eq.([21) and 
fi ~ (7^ + 1. We assume that here as everywhere throughout the paper the density fluctuations (the total 
one and the difference between the layers) are incommensurate with the lattice. In the context of model 
([3]) it is achieved by a proper choice of the chemical potential. 



Following Schulz [l3| we will treat this model at low energies using bosonization technique (see als6 
Ref.fl^). The model describing each chain is the spin S=l/2 XXZ model; in the continuous limit it is 
equivalent to the Gaussian model: 



(4) 



where a = 1, 2 labels the tubes (layers) and Qa is the field dual to [dx'da{x),<^a{y)\ = -'in6{x — y). 
The Luttinger parameter K is determined by the intra-chain interactions, that is, in the case of the model 
([3]) by the ratio Vo/t|[. We assume that > 0; then the continuum limit of the spin operators is given 
by the following bosonization formulae: 



lV27r*„ 



(2^ao)V2 



C 



gi\/27r(ea+*a)+2iA;FX _|_ gi\/27r(-eo+$a)-2iA;F2; 



y/n (27rao)^/^ 



(5) 



where dots stand for less relevant operators and C, are amplitudes determined by the short range 
physics and qq is the short range cut-off. The Fermi momentum kp for each chain is determined by its 
chemical potential fi with the convention that kp(fi = 0) = 0. As we have mentioned above, we will 
always assume that fi ^ so that the spin fluctuations are incommensurate with the lattice. 
Substituting ([5]) into ([3|) and defining the fields 



$1,2 

ei,2 



[dxQaix), $fc(y)] = -iSab5{x - y), 
we obtain the Hamiltonian H = H-^- + iJ_ , where Hj^ describes the symmetric mode 



Ax 



and i7_ contains only the anti-symmetric fields: 



H-= da;|^[(a,$_)2 + (a,e_)2] + COS (yi;^:^*^) -l/iCos(V4^/if-e_)|, 



(6) 



(7) 



(8) 



where Vi ~ Vi , t ~ t and 



K+=K± 



Vi 

2lTV 



(9) 



At this juncture we note that the Hamiltonian ([7]) describes a mode which is always critical. The 
corresponding order parameter is = crj^'cr^; according to ([5]) 



(10) 



Although in one dimension such order parameters with continuous symmetry do not have vacuum aver- 
ages, at T = their correlation functions exhibit slow (algebraic) decay. In the present case 



< ^'+(x, 0)^'^(.t', 0) >- \x - X 



'1-2X4. 



(11) 



However, there may be operators with correlators decaying faster than that of 4*+. They are different in 
different phases of our model. Depending on which of the cosines in ([S]) takes over, the ground state of 



this model describes either quasi long range superfluid order in each tube or pair density wave state. The 
latter state has a singularity in the density-density correlation function at the finite wave vector 2kp. 
When both cosines are relevant (that is at 1/2 < K- < 2) these states are separated by a quantum critical 
point (QCP), the location of which is approximately determined by the relation (Vj/A)^- ^ {Vc/Ay/^- , 
where A is the ultraviolet cut-off. Fulfillment of the above condition on if _ is essential for the subsequent 
arguments. 

The vicinity of the QCP can be studied analytically when w 1 (this will be our assumption 
throughout the rest of the paper). In that case it is convenient to refermionize with the result 



Ativ-{K^ - l)pRPL'ilRilL + 2im+pRPL + 2im^r]R-qL^, 



(12) 



Where ■m± 
fermions: 



t±^ ± Vi and pl,r. and ?7L,fl are left- and right-moving components of Majorana (real) 

1 



PR,l = 



VR,L = 



1 



: COS 



v^($„±e_) 
V^($_ ±e_) 



(13) 



We comment that the fcrmionization of the system in terms of real (Majorana) fermions is consistent 
with the Ising type symmetry breaking. This model is equivalent to the continuum limit of two quantum 
Ising (QI) models coupled by the energy density operators [l^. The transition occurs when one of the 
Majorana masses becomes zero. To access the correlation functions we need to express the original spin 
operators in terms of the Ising model fields: 



/2^ 



-.dxQ+ + c 



.(.14) 



where Q ~ tt/oo -I- 2kp, s± {p±) arc Ising order (disorder) parameters for the Ising models represented 
by p and 77 fermions respectively. These operators are nonlocal in terms of fermioms. Their explicit 
expressions are not needed here, all we need to know is that in the part of the phase diagram m > we 
have (a) 0, (p) = and for m < we have (a) = 0, (p) ^ 0. Therefore at tj^ > Vi > when both 
masses have the same sign both a± have vacuum averages. Replacing these operators in (jl4p by their 
vacuum averages we get 



(15) 



Since has a gapless spectrum, the corresponding correlation function decays algebraically with in the 
exponent K+/2 which is four times smaller than the exponent for 5'+ (fTO|) . In the other phase m+m- < 
and a similar replacement can be done for the oscillatory part of the density operator yielding 



al ^ ±cos(Qx -I- y/T:Q+)[< s >+< p >_ 



so that 



<<(x,0)a,^(x',0) >r 



cos[Q(a; — x')] 



(16) 



(17) 



The latter situation corresponds to Pair Density Wave (PDW). The density oscillations in this phase exist 
alongside the superfluidity of coupled pairs, though in the presence of disorder PDW is pinned [iBl and 
the superfluidity is not destroyed. The QCP separating the two phases occurs when one of the masses (if 
t± > 0,Vi > it is always 77i_) becomes zero. Thus, for two chains N = 2 the transition occurs in one 
Ising model and since Matsubara time correlation functions of quantum Ising model at T = are the 
same as the correlation functions of 2D classical Ising model, it belongs to the universality class of d = 2 
classical Ising model. The correlation length exponent \s v = 1. 



B. N tubes symmetrically coupled 
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Now we consider a system of iV > 2 tubes coupled to each other in such a way that each tube 
interacts with aU others. The treatment in this case is more comphcatcd since we have to resort to 
non-Abehan bosonization (see, for instance, [l3i[I3)- ^ starting point we take non- interacting tubes 
with SU(2) symmetry (Vb = t\\). Then an invidual tube is equivalent to spin S=l/2 isotropic Heisenberg 
antiferromagnet and in the continuum limit is described by the SUi(2) Wess-Zumino-Novikov-Witten 
(WZNW) model. The sum of TV SUi (N) WZNW Hamiltonians can be decomposed as (see, for instance 

m 

SUi{2) + ...SUi{2) = \SU2{N)/U'^^\l)] X Znx U{1) (18) 

This decomposition should be understood in the sense that operators (primary fields) of the critical theory 
on the left hand side of the identity can be written as products of operators belonging to the critical field 
theories on the right hand side. Decompositions of that kind can be very helful outside criticality if the 
perturbation happens to be such that it does not act in at least one of the sectors. Decomposition ^TE\i 
is consistent with the fact that central charges of the theories on the left- and right hand side of (|18p are 
equal: 



N : 



2(7V2 - 1) 
N + 2 



(A^-l)] + (^-l)+l. (19) 



The U(l) subsector of ([T5)) corresponds to the symmetric bosonic phase (N-tube generalization of $+ 
from the previous subsection). Since the inter-tube interaction does not contain this field, it remains 
gaplcss. As far as the other sectors are concerned, we will leave a detailed analysis to future publications 
and only formulate some conjectures. Information extracted from our numerical calculations suggests the 

following scenario. Close to the critical point the SU2{N)/U^^^{1) -sector (Gepncr's parafcrmions [l9| ') 

is weakly coupled to the rest. This coset sector remains massive throughout the entire phase diagram, 
at least in the part where t± > 0, Vi > 0. Its analog for iV = 2 is p Majorana fermion. The Zjv sector 
is the one where the critical point is located. The relevant perturbation around the critical point can 
be guessed from the numerics which yields i/ = 5/6 ~ 0.833 for = 3 and ly — 0.75 for A^ = 4. Using 
the relation v = 1/(2 — d), where d is scaling dimension of the operator responsible for the deviation 
from criticality, we find d « 0.8 for A^ = 3 and d « 0.67 for A^ = 4. On the other hand in the model 
of Zjv parafermions there is an operator with scaling dimension d = 4/ (A^ + 2)[l^ [2^ which reproduces 
perfectly the numerical values of d (see below) . 

IV. NUMERICAL RESULTS 

The QPT transition discussed above is characterized by disappearance of the ODLRO in all M-body 
density matrices, where M = 1,2, ...^ N — 1. Specifically, as the inter-layer interaction is increasing the 
ODLRO existing in all M-body density matrices must eventually vanish up to the order M ^ N — 1. At 
the same time, ODLRO remains, practically, unaffected in the N-body density matrix. 

A. M-body density matrix 

The M-body density matrix Dm can be written explicitly as 

Dm{{{xi,zi),...,{xm,zm)};{{x'i,z[),...,{x'm,Zm)}) = { ]J ai„,zm IT "^™'^,r.')' (20) 

m=l,...,M m' = l,...,M 

where (...) stands for the quantum-thermal averaging. 
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FIG. 2: (Color online) d{G'j,)/dU versus the interaction strength [/ for L = 60,100,200,300 with /3 = i. Inset: 
{G'i) versus U for L — 100. The SF phase corresponds to (G3) ^ 1 and the CSF to (G3) « 0. The transition 
point SF-CSF for a given size can be identified by the maximum of dGz/dU reaching the thermodynamics limit 
at Uc ~ 0.61. 



In ID SF Dii ') ^ l/\x — x'l'', b < 1, exhibits algebraic order at large \x — x'\. In the CSF, 

Di{x, z;x' , z') ^ exp{—\x — x'\/£,o), ^0 ~ Ij that is, it becomes short ranged at T = regardless of the 
filling factor. Thus, the TV-body density matrix, on one hand, is characterized by the exponential decay 
Dn{xi, Xm.]x'i, x'jy^) ^ exp(— Ixmj — cCmj |/Co) wlth respect to any pair Xm^ , Xm^ of coordinates from 
the set xi, {or x'l, ...,x'^). On the other hand, there is the algebraic order £) at ^ 1 / \Rcm~ R'^^\'^ , c > 
0, with respect to the "center of mass" coordinates Rem = [xi + ... + xn]/N and R'„n = [x[ + ... + x']^]/N 
defined, respectively, for the first and the second sets of the coordinates as introduced in Ea. ipO)) . provided 
\Rcrn - Xm\ < ^0 and - x'^\ < for aU to. 

The transition SF to CSF can be detected by critical behavior of any density matrix. In particular, the 
same criticality controls the long-distance behavior of Dj\j with respect to \Rcm ~ Xm\ (or |i?cm ~ ^ml)- 
That is, in SF phase is trivially long-ranged with respect to \Rcm ~ Xm\ because can simply be 
factorized into a product of Di. In contrast, in the CSF-phase, while exhibiting ODLRO with respect to 
Rem — R'cm^ IS short-raugcd with respect to \Rcm — Xm\ (or |i?cm ~ ^ml)- Thus, the criticality can 
also be detected by measuring the behavior of the relative distances Xm (or x'^). Specifically, we have 
considered the square of so called gyration radius [lo| as the mean of 



i?2 



^ ^ ~ ^"]^ (21) 

m,n=l,2,...,N 

with respect to the first set of the coordinates of 13 at defined in Eq. ([20|) . provided the coordinates from the 



■9 ^2 

m,n=l,2,...,N 

yi Dn define( 

second set are kept within some distance ^ from i?^^ 2l[. More specifically, Xk, where fc = 1, 2, .., A^, 
represents the x-coordinate in the A;— th tube. 

^0 ~ 4(3-6) ' 

In what follows we will be calculating the mean of the ratio Gjv = Rg/R-o, so that it is changing from 
Gn ~ 1 in the SF state to Gn ^ 1/L^ in the CSF phase. It is worth mentioning that Rg can be 
viewed as a typical width of a chain. For strongly bound case this width is ^ ~ Ij arid in the SF phase 
it is ^ L, and, thus, it exhibits critical behavior typical for correlation length. 



In the SF of a length L, R^ ^ R^ ^ IT^L^ ~ 0{L^), and in the CSF RI^f ^ C'(l) ~ Co << 



B. J-current formulation 



Hamiltonian ([2]) and the gyration radius (j21|) have been used for ab initio simulations of a single chain 
(with exactly one polar particle per layer) for the case t± = It was found that the chain can 
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0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 <G3> 

FIG. 3: (Color online) d{Gci)/dU versus (G3) for sizes L — 20, 40, 300 rescaled by a factor A(L) in order to 
achieve collapse to the curve L — 100 (A(IOO) = 1). Inset: d{Gs)/dU versus (G3) for the same sizes. 



v=0,835+-0,015 



32 64 128 256 



FIG. 4: (Color online) The rescaling factor A ^(L) versus L for A'' = 3 from Fig[3] The slope gives the correlation 
length exponent v = 0.835 ± 0.015. 



P(E) 



0,80 0,82 0,84 0,86 0,88 E 



FIG. 5: (Color online) Bimodal histogram of energy P{E) for TV = 5 tubes with L = 400, 13 = 400. The first-order 
transition happens at Uc ~ 0.7235. 
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6-5 Ln L 



FIG. 6: (Color online) The maximum value of d{GxY)/dU versus InL for the case of zero interlayer tunneling, 
A'^ = 4, K{z) — oo. The solid red line is the fit by the finite size scaling ansatz for the BKT-transition: 
d{GxY)/dU = A\n^{L/Lo), A = 0.205, L„ = 1.98. 



undergo quantum roughening transition with the tuning parameter being the interaction strength Vd- 
The transition is, practically, insensitive to the interaction range. 

The Monte Carlo simulations at finite densities n < 1 ( incommensurate with the lattice along the 
tubes) in each layer have been conducted in the discrete-time J-current-typc formulation of the 
Hamiltonian ([2]) [Ifll • For the purpose of analyzing the universality of the transition this approach turns 
out to be much more efficient than the ab initio one. Here we will be using the model where the inter- 
layer tunneling is allowed. The actual dipolc-dipole interaction will be replaced by onsite attraction 
between neighboring layers, with the intra-layer dipole-dipole repulsion ignored. This approximation 
becomes essentially exact when n << 1: while the inter-layer attraction is not affected, the intra-layer 
dipole-dipole repulsion between atoms scales as '--^ n"^ — > in Id and, thus, becomes irrelevant. The 
corresponding space-time action, then, becomes 



(22) 



where Jh is the integer bond current obeying Kirchhoff's conservation law |22| . In some sense, these 
conserved currents represent world-lines of particles in imaginary discrete space-time, with Hj being the 
action in Feynman's path integral. The summation in (|22p is performed over all space-time bonds b 
(coming out from a space-time site (x, r, z) either along ±x or along imaginary time ±f or along ±z 
directions); V^Jj, = Jb{x,T,z -|- 1) — Jb{x,T, z); fi denotes chemical potential. [Here we tuned fi to have 
1/2 filling of bosons per site in each tube]. Periodic boundary conditions along space 0<a;<L— 1,0< 
z < N — 1 and along imaginary time < r < /3, with (3 = L, where L = 2,3, have been used. The 
coefficients K,U can be related to tz,z' and Vxz-x'z' from Eq.(I2]): K{z) sa l/t_L, K{x) = K{t) ss 
U w Vd- The case K{z) = oo corresponds to zero inter-layer tunneling (studied in Ref . jTo| ) . Here we 
will focus on K{z) — K{x) — K{t) situation as the one which naturally represents the whole universality 
class. 

It is worth emphasizing that the inter-layer attraction between two particles (in neigboring tubes) 
located, respectively, at the space-time points {x,z,t) and {x,z ± 1,t) is described by the terms ~ 
~U J{x, z,t)J{x, z ± 1,t). Accordingly, when these two particles form a bound state, their world-lines 
stay close to each other to gain the binding energy ~ U . Similarly, a chain of several particles is 
represented by a bundle of several world-lines forming a membrane in the space-time. 

The action can be viewed as a coarse grained dual representation of the Hamiltonian While 
being not precise for quantifying finite energy (non-universal) properties of the system, the J-current 
model [1^ belongs to the same universality class as the original model Thus, for the purpose of this 
work and for sake of numerical practicality, it will be sufficient to study the model 



We also note that, while being formally defined on the lattice, the model © and its dual formulatiifi 
account well for the continuous space-time situations at low energies as long as the filling factor n 
remains incommensurate with the lattice. Long-range intra-layer repulsion may complicate the situation 
by inducing crystalization at, say, n ~ 0.5 and, thus, shifting the CSF phase to lower densities. Such 
feature, however, does not affect the universality of the N-SF to CSF transition, and, in order to establish 
it in a most efficient way wc simply turn off the intra-laycr repulsion and study the case n = 0.5. 

Monte-Carlo simulations of the model ((22|) have been performed within the Worm Algorithm ap- 
proach [23| . Green's function in imaginary time (as well as the density matrix (j20|) ') is given by the 
statistics Dpf({xm,T„i, Zm}; {a^m, t^,, z'j^^}) of "sources" and "sinks" of the bond currents located, respec- 
tively, at (xm,T„i, z,n), m = 1,2,...,N, and (x^,t^j,z^), m = l,2,...,iV, lattice points. In order to 
insure the condition li?^™ — a^ml ^ ^Oj while {xm,Tm, Zm), m = 1, 2, N, are free to take any value, we 

have convoluted DN{{Xm,Tm, Zm}; {x'^,T'^, Z'^}) with P = CXp(- Xlm.nfl^^m - a^'nl + k'm ~ <l]/Co) aS 

DN{{xm,Tm, Zjn}; R'cm) = / Dx' Dt' D z' D N PS (R'^^ ~ J2m Kn/^) J ^ud, accordingly have evaluated the 
means of the normalized gyration radius (Gjv) and of the center of mass distance {\Rcm — R'cml) where 
(...) = Z-^ J DxDTDzdR'„^...DN, Z = J DxDtDzcLR'^^Dn ■ 

For sake of numerical efficiency we have symmetrized the model \22\ by choosing U{h) independent 
of the type of a bond, that is, U{h) = U . The CSF phase has been identified by the condition {\Rcm — 
R'cm\)/L = const and (Gtv) ~ o{L~^) for U > Uc, where Uc corresponds to the QCP. In the SF phase (that 
is, U < Uc), while the first condition remained, practically, unchanged, (Gjv) ~ 1 with high accuracy. 
The criticality of the SF-CSF transition has been analyzed through evaluating the divergent behavior of 
d{GN)/dU in the vicinity oi U = Uc- 

C. Finite size scaling of the gyration radius 

As discussed above, d{GN)/dU exhibits singularity in the limit L oo. The change from {Gn) ~ 1 
to {Gn) ~ occurs in a narrow range 5U = \U — Uc\ around the QCP Uc ^ 1- Such a behavior is clearly 
seen in Fig. [2] the range of the transition 5U narrows as L increases. 

This range is controlled by the diverging correlation length ^(C/) ~ \U — Uc\^'' , v > Q. According to the 
finite size scahng approach, {Gjq) can be represented as some regular function F{y), y — L/£^{U) varying 
from F{y = 0) = 1 to F{y = oo) = over the range y - 1. Thus, d{GN)/dU w F'y/SU - L^/''. Loosely 
speaking, one can view this relation as d{GN)/dU ~ 1/6U, 6U w L~^/^ — >• 0. 

We have evaluated this derivative numerically by Monte Carlo [2^ and constructed the graphs 
d{GM)/dU versus (Gat) by scanning over U around the critical point Uc for sizes L = (3 = 10,20, ...300. 
These graphs turn out to be self-similar so that d{GN) /dU for all sizes > 10 collapsed on a single master 
curve by simple rescaling of d{GN)/dU for size Li to another size L2 as d{GN)/dU — \{L)d{G m) / dU . 
Then, the rescaling coefficient A, which represents the inverse width 5U as A oc 1/5U , has been plotted 
in the log-logL axes in order to determine the critical exponent v. The results of this procedure are 
presented on Figs. I3I4I for the case iV = 3. The same procedure has been used in the cases TV = 2, 4 as 
well. The found exponents are: v = 0.972±0.02 for TV = 2, ^ = 0.835±0.015 for = 3, = 0.735±0.015 
for A = 4. [The errors include statistical errors as well as the systematic errors due to the subleading 
contributions] . We note that the value of v for A^ = 2 is consistent with the d = 2 Ising (or q = 2 Potts) 
universality. We also note that the values of for A^ = 3 and A^ = 4 are consistent with the corresponding 
ones V = 0.837 and v = 0.756 obtained by the Renormalization Group calculations for the g = 3,4 2d 
Potts model [13 . 

For A^ > 4, the transition was found to be of first order. It has been detected by the bimodality of 
energy histogram, Fig. [5j While for A = 5 such bimodality develops on sizes L > 400, /3 = L, for A = 8 
it is already well developed at L = 160, f3 = L. 

The finite size analysis has been applied to the case of zero inter-layer tunneling as well, when the 
transition is expected to be in the BKT universality. That is, ^ ~ exp(...|J7 — Uc\^^^^). The variation of 
the gyration radius (Gxy) in this case can also be represented by some regular function F{y) characterized 
by the range y - 1 with y = L/C Thus, d{GxY) /dU w F'y - (SU)'^^^ - (ln(L/Lo))3 (where stands 
for some microscopic scale) at its maximum. The maximum value of this derivative has been plotted 



as a function of L = 30, 600 in FiglS] As can be seen the fit of the data is consistent with the \v?^i: 
dependence with high accuracy. 



V. SUMMARY AND DISCUSSION 



As shown above, the Majorana fcrmions description of the Ising-type transition to CSF state in the two- 
chain system is consistent with the numerical evaluation of the correlation length exponent v — 0.972±0.02 
(versus its exact value v = 1). The quantum transitions in the cases A'^ = 3, 4 are characterized by enlarged 
emerging symmetries: Z3 for A = 3 and Z4 for A = 4. The predicted values = 5/6 « 0.833, A^ = 3, 
and V = i/ A = 0.75, A^ = 4, arc matched well by the corresponding numerical ones v = 0.835 ± 0.015 
for A^ = 3 and v = 0.735 ± 0.015 for A^ = 4. The symmetry enlargement occurs despite the short-range 
nature of the tunneling between the layers. In other words, the critical behavior proceeds as though the 
tunneling between all layers is the same. Such feature — identical interaction between all elements — is 
typical for the standard Potts model [ll| and should be contrasted with the p-clock model. 

The detection of the CSF order as well as the criticality to SF state can be based on measuring 
field-correlators. In the case of the A'-layered structures supporting indirect excitons [H] this means 
analyzing the excitonic emission of light. Indeed, as we mentioned in the Introduction, such emission 
is electric-dipole active. Thus, the electric field from zth layer in the far zone at point f is E{r) ~ 
J dx ex.'p{ikxx)ili{z , x) , where k is the wavevector of the light emitted in a given direction k/k from the 
excitonic sample to the point f and ip{z,x) stands for the excitonic field. Thus, in the CSF phase where 
{ijj^{z,x)4'{z,0)) is short ranged (with respect to x), the emission of light will be essentially the same 
as from a thermal source, that is, incoherent and isotropic. Conversely, emergence of the SF order 
will be characterized by narrowing of the emission. Meanwhile the amplitude of emitting N-photons is 
determined by the relation ~ E(ri)...E{rM) ^ tp{l,x)...tp{N,x) ~ ^n, where 'i']^, Eq.(IT]), carries the 
information about the CSF phase. Thus, the A^-photon emission correlated measurements can reveal 
the CSF order. We will consider specific proposals for detecting the CSF order as well as the discussed 
criticality in greater detail elsewhere. 

When this work was prepared for publication we learned about the preprint by Lecheminant and Nonne 
(25I which results have a substantial overlap with ours. 
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